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ABSTRACT 

Low-energy cosmic rays (CRs) are confined by self-generated MHD waves in the mostly neutral 
ISM. We show that the CR transport equation can be expressed as a continuity equation for the CR 
number density involving an effective convection velocity. Assuming balance between wave growth 

and ion-neutral damping, this equation gives a steady-state condition n CI oc n ; up to a critical 
density for free streaming. This relation naturally accounts for the heretofore unexplained difference 
in CR ionization rates derived for dense diffuse clouds (McCall et al. 2003) and dark clouds, and 
predicts large spatial variations in the CR heating rate and pressure. 


Subject headings: ISM: cosmic rays 
1. INTRODUCTION 

Low-energy cosmic rays (CRs) are coupled to a num¬ 
ber of important astrophysical processes. Most of the 
integrated interstellar CR number density, energy den¬ 
sity, pressure and ionization rate are contributed by the 
low-energy part of the spectrum because of its steep¬ 
ness (see compilation in Antoni et al. 2004). 1 10 GeV 
CRs control the ionization fraction in the Earth’s lower 
atmosphere, and so cloud formation and lightning pro¬ 
duction (e.g. IStozhkov) 12003 1. Tropospheric and strato¬ 
spheric chemistry are affected, ev en at energies as low 
as ~ 10 MeV (e.g. ICrutzen et alJi|1975 l. In protostellar 
disks, ionization by ambient low-energy CRs may reg- 
ul ate dis k tur bulence and affect planet formation (e.g. 
iMatsumura fc Pudrhdf2005 1. On the scale of the inter¬ 
stellar medium, low-energy CRs of ~100 MeV (Webber 
1998) dominate the ionization fraction and heating of 
cool neut ral gas, especia lly dark UV -shielded molecular 
regions llColdsmith fc Langer 111978th On larger scales, 
CR “pressure”, most of which is contributed by low- 
energy CRs, may help confi ne the Galactic disk an d drive 
the Parker instability Isee lHanasz fc Leschll200 0l. which 
may itself contribute to the formation of large condensa- 
tions that for m stars, and even drive a Galactic dynamo 
(lParkeml992lk t he CR pressure may als o affect the hot 
coronal ISM ( ISchlickeiser fc Lerchdl19851) . For these rea¬ 
sons, any strong spatial variations of the CR number 
density could lead to important thermal, chemical, and 
dynamical effects. 

In almost all previous work it is assumed that the CR 
density and spectrum do not vary significantly in space 
for length scales smaller than the scale of variation in 
space dens ity of CR sources or ionization and spa llation 
losses (e.g. iHunter et al .11 19971: iWolfire et al .112003th The 
primary rationale is that a superposition of CRs prop¬ 
agating diffusively to a point in the Galaxy from many 
stochastic sources (e.g. supernova remnants) gives an 
rms variation in CR density of order 1% for typical values 
of parameters (Lee 1979; see Berezinskii et al. 1990, sec. 
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III. 10). A more detailed calculation including the spa¬ 
tial distribution and discreteness of sources gives varia¬ 
tions that are typically less than 10-20% (Busching et al. 
2005). This argument is not valid for low energy CRs be¬ 
cause it neglects the possibility that self-confinement of 
CRs (see below) or ionization losses greatly reduce their 
propagation distance from the sources. The evidence for 
CR h omogeneity inferred from EGRET 7-ray data (e.g. 
IDigel et aHl200lf) applies only to a spatial resolution of 
approximately 5 degrees and to CR energies larger than 
of interest here. Furthermore, the 7-ray emissivity is de¬ 
rived from an integral over long lines of sight, to which 
individual clouds or cores may be smal l contr ibuti ons, e s- 
pecially at low Galactic latitudes (see IAharonianH[200lll . 

The purpose of this Letter is to show that large spatial 
variations of the CR number density should indeed exist 
in the ISM, based on the standard CR transport equation 
(§ 2), if low-energy CRs are confined along flux tubes 
by scattering from self- generated Alfvcn waves (§ 3). 
ISkilling fc Strand (1197fif) examined a model to exclude 
CRs from molecular clouds based on screening the CRs 
by ionization losses enhanced by CR self-confinement, 
but their model assumes that CR self-confinement sets 
in suddenly at the edges of molecular clouds because of 
an extremely large assumed cloud column density. The 
process examined here is completely independent of that 
model. Multiple magnetic mirrors could also lead to CR 
variations of a factor of a few , depending on adopted pa¬ 
rameters dCesarskv fc Volld 1 978|) , a process that could 
be more important in the presence of tangled fields. In 
reality several effects may contribute to variations on 
roughly the same scale, but the effect found here has 
not been previously recognized, and is capable of giv¬ 
ing CR variations up to two orders of magnitude. We 
show that, in a steady state, the CR density should lo¬ 
cally scale with the square root of the ion density, up 
to densities above which damping of the waves allows 
the CRs to stream freely (§ 3), giving a sharp decline 
at densities around 500 cm -3 (for 100 MeV protons), 
typical of the transition from diffuse gas to dark molec¬ 
ular clouds. These variations may explain the apparent 
discrepanc y b etween the large ionization rate derived by 
iMcCall et alJ ( 2003 1 for the diffuse region along the line 
of sight to C Persei and the ionization rate estimated in 
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dark molecular clouds (§ 4). They should also result in 
large variations in CR pressure and heating rate in the 
ISM and elsewhere. 


2. CONTINUITY EQUATION FOR SELF-CONFINED 
COSMIC RAYS 


It is well -known that CRs streaming along field lines 
at a s peed larger than the Alfven speed generate MHD 
waves llLeichell 19671 I Wentzelll 196811 . The CRs in turn in¬ 
teract with these self generated waves through resonant 
pitch-angle scattering with waves whose wavenumbers 
are multiples of t he particle gy rora dius for a particle of a 
giv en ener gy (seelCesarskv & Kulsrudlll973D . As shown 
bv ISkillind ( 1971 1 and others, the waves keep the CRs 
confined to stream at a speed only slightly larger than 
the Alfven speed, if the waves are not efficiently damped. 
Models of Galactic CR propagation that rely on this ef¬ 
fect in order to increase the escape time of CRs from 
the Galaxy are called self- confi nement models, and are 
reviewed in iWentzell (|1974f> and iCesars EUIlU. Self¬ 
confinement is generally believed to be efficient for CRs 
of energies less than about 1 00 GeV in mostly neutral 
interstellar material iKulsrud fc CesarskvHl97T ). 

To demonstrate how self-confinement can lead to large 
CR density variations, we begin with the CR transport 
equation. The usual kinetic equation for the phase space 
distribution, /(x,p, t), of CRs interacting with a back¬ 
ground plasma or self-generated MHD wave field can 
be derived from the Vlasov equation, coupled with the 
CR particle equation of motion and Maxwell’s equations. 
A number of reasonable assumptions allows transforma¬ 
tion of this equation into an equation for the distribution 
function of guiding centers, assuming the CR anisotropy 
is small. This equation was derived and d iscussed i n var¬ 
ious fo rms by iKulsrud k Pearcd (I1969H: ISkillind (F - 


ious torms d v iivuisrua__ _ _ _ 

[HzpLSlIm and many others. iLu V Za.nH ll2flfH 
and lLu et alJ (| 2002 |f give a useful collection of references 
(see Berezinskii et al. 1990, Schlickeiser 2002 for detailed 
derivations). 

Neglecting diffusion perpendicular to the magnetic 
field, introducing a non-relativistic background medium 
velocity with a Galileian transformation (V <C c) and 
including continuous momentum losses, the result is 

f + v.v/-v.[Mx,p)v/] + |v.v| = 


J _d_ 

p 2 dp 




,dp 
d t' 


S(x,p,t) 


( 1 ) 


This well-known transport equation, neglecting the 
terms on the right, was derived phenomenologically by 
I Parker : 196 ') 1. It is especiall y com mon in studies of helio - 
spheric cosmic ray transport llFerreira & Potgieteil2004lf . 

On the lhs the second term is the convection of / by 
the background plasma bulk motion at velocity V, which 
can include guiding center drift velocities (important in 
heliospheric transport but not in most ISM conditions). 
When this background motion is due to MHD waves, the 
appropriate velocity depends on the distribution of the 
directions of wave propagation relative to the streaming 
of the CRs. If the waves are self-generated, as we as¬ 
sume here, the directions of the waves and CRs are the 
same, and V can be re placed by the Alfven speed, Va 
llBerezinskii et al1 ll990. ch. 10). The third term repre¬ 
sents the interaction of the CRs with the MHD waves 


in the diffusion approximation; k(x,p) is the pitch-angle 
averaged spatial diffusion coefficient, which derives from 
the slight anisotropy of the CR distribution. The fourth 
term, often called the adiabatic or the Compton-Getting 
term, represent momentum convection. On the rhs the 
first term represents momentum (or energy) diffusion, 
the second term continuous momentum losses due to in¬ 
teractions with plasma particles (e.g. spallation, ioniza¬ 
tion and radiation losses), and S(x,p,t) is the source 
distribution function. 

We can safely neglect all the terms on the rhs. It is 
well known that momentum or energy diffusion is slow 
compared to pitch-angle diffusion (transformed into the 
spatial diffusion on the lhs), by a factor of order Va/c 
(see Berezinskii et al. 1990, ch. 10). We neglect ioniza¬ 
tion losses because they require a column density of order 
100 g/cm 2 for 100 MeV protons, corresponding to large 
length scales even when tangled fields or self-confinement 
are taken into account. The sources (e.g. supernova 
remnants) can be assumed far from the relatively small 
(^0.01 to 10 pc) volume under consideration. Similarly, 
we assume the mean magnetic field varies only over scales 
much larger than the scattering mean free path, so that 
we neglect mirroring and drifts due to mean field varia¬ 
tions (weak focusing limit). Our resulting scale of varia¬ 
tions is probably comparable to that of field variations, 
so a full calculation should include this focusing term. 

Neglecting all terms on the rhs, rearranging the advec- 
tion term with the Compton-Getting term, and multi¬ 
plying by Airp 2 , equation 0 yelds 

If + V • (V 5 ) - V • [k(x,p)Vg] = iv • (pg) (2) 

where g(x, p , t ) = 47rp 2 /(x, p 1 t) is the differential number 
density of CR particles. An integration over p then gives 
(assuming that p g vanishes at infinity) 

r)n — 

+ V ■ [V n cr - fc(x)Vn cr ] = 0 (3) 

where n CI = / 0 °° g(x,p,t)dp is the total num¬ 

ber density of CR particles, and k is the 
momentum-averaged spatial diffusion coefficient, 
fc ( x ) = r k(x,p)g(x,p,t)dp/n cr (x,t). 

Equation m can be cast in the form of a continuity 
equation by noticing that the diffusion along the field 
is the divergence of a flux that can be represented by 
the product of a diffusive streaming speed, Vdig, and the 
CR number density, n cr . We can then define an effective 
CR streaming speed, v st , as the sum of the convection 
velocity, V, and the diffusive streaming speed, Vdiff. In 
that case, an equation expressing a steady state for the 
CR number density in an Eulerian frame is 


V • (v st n CI ) = 0 (4) 

A similar connection between CR transport and a con¬ 
tinuity equation has been noted, in different contexts, 
by Skilling (1971, 1975a), Earl (1974), Schlickeiser & 
Lerche (1985), Beeck & Wibberenz (1986) and Bieber 
(1987). Integrating equation Q) over a small volume, 
using the divergence theorem with a closed surface that 
corresponds to a segment of a magnetic flux tube, and 
using the inverse proportionality between the flux tube 
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cross section and the magnetic field strength, B , we ob¬ 
tain 


i' s t n CT 

~B 


= const 


( 5 ) 


As significant spatial variations of both B and v s t are 
certainly present in the ISM, we should also expect sig¬ 
nificant variations in the CR number density, n CI , CR 
pressure, P cr , and CR ionization rate, £ cr , as discussed 
in § 4. 


3 . COSMIC RAY STREAMING VELOCITY 

We can show how n CI and associated quantities should 
vary with ISM parameters by deriving an expression 
for t) s t. If the CR scattering is primarily due to res¬ 
onant scattering off magnetic waves generated by the 
CRs themselves, v s t can be computed by requiring that 
the wave growth rate is ba lanced by the wave damp¬ 
ing rate (see ilWentzeli 1974 . sec. 2.3-2.5). Consider¬ 
ing only protons, the growth rate of waves propagat¬ 
ing in the direction of the magnetic field, B, as a func¬ 
tion of the CR stre aming velocity along B, v 8 t is (eg 
iKulsrud fc Cesarskvlll971f) : 


_ tt (7 - 3)fl 0 n>(k z ) rn H f 3 _ \ 

4(7-2) m mi \7 Va J 


( 6 ) 


where flo = eB/mnc is the non-relativistic cyclotron 
frequency of a proton of mass mn, the CR spectrum 
has been assumed to be a power law in momentum with 
exponent 7 (empirically 7 = 4.7 for E > 10 GeV), mi 
is the ion mass, and n> r (h z ) is the number density of 
protons with momentum p > eB/k z c, corresponding to 
the resonant condition. 

In the mostly neutral ISM damping is primar- 
ily due to collisions o f ions with neutral particles 
iKulsrud fc Pearc3 '!i 969(l . because neutrals do not take 
part in the wave motion, as the wave frequency is larger 
than the ion-neutral collision frequency at the scales of 
interest . 3 The ion-neutral damping rate is then: 

„ 1 m „ 

Tin = -n n (0V n ) in — ( 7 ) 

z rri[ 

We take the collision rate (av)i n = 2.1 x 10 -9 cm 3 s _1 in 
diffuse regions and (crv)^ = 1 . 6 xl 0 -9 cm 3 s _1 in molec¬ 
ular regions llOsterbrockllQfilH . assuming n(He)/[n(H)+ 
2n(H 2 )] = 0.14. 

Assuming a balance of wave growth and damping, 
T = Tin (such balance is reached in less than about a 
year for typical parameters in the absence of nonlinear 
wave cascades; see Kulsrud & Pearce 1969), equations © 
and Q {j ive an expression for the CR streaming velocity 


v s t {p) 


iVa 2 (7 - 2 ) n n {av n ) m m n 
3 [ n (7 - 3) fi 0 n>(p) m H 


( 8 ) 


The importance of the second term on the r.h.s. of equa¬ 
tion @ depends on the gas density and fractional ioniza¬ 
tion. For mostly neutral ISM conditions and for ionizing 


3 Nonlinear cascade damping of hyd romagnetic waves (e.g. 
ISkillingll975bl : lFarmer fc GoldreichM2Q04D is slow compared to col- 
lisional damping in the cool mostly-neutral ISM. 



Fig. 1.— CR density versus gas density predicted for mostly 
neutral diffuse gas and dark cloud cores by iteratively solving equa¬ 
tions 0 and Jo). Details are given in the text. 


protons of energy E < 100 MeV, the second term is <C 1 
up to a gas density nn.fs ~ 500(Bo/10 pG) cm -3 , and 
low energy CRs are well confined, streaming at a veloc¬ 
ity of order the Alfven speed. In this case, equation © 
yields: 


1/2 

n cr oc n ; 


(9) 


Notice that the magnetic field strength has dropped 
out for this regime; it does affect the value of nn,fs above 
which eq. © no longer holds. As the gas density is 
increased above nn,fsj the second term on the r.h.s. of 
equation © becomes important and the CR streaming 
velocity approaches the particle velocity, causing a drop 
in the value of n cr , based on equation ©• 


4 . RESULTS AND CONCLUSIONS 

Figure 1 shows the result, for 100 MeV protons, of 
iteratively solving equations © and © for the depen¬ 
dence of total CR number density on total gas density, 
takin g n> = 7 y r . W e have assumed «; = 1.4 x 10 _4 n 
llCardelli et all I1996H and B = B 0 in diffuse clouds, 
and m = 10 _ 4 (n/ 10 4 cm“ 3 ) 1 / 2 (C C r /2 x lO^V 1 ) 1 / 2 

Bo(n /200 cm -3 ) 1 / 2 
in molecular clouds, 
with B 0 = 10 pG. These are crude approximations due 
to the large observed scatter. The vertical normalization 
assumes a CR energy density of approximately 1 eV/cm 3 
at n = 0.1 cm -3 , correspondin g to dem odulation of the 
energy spectrum near the Sun llWebbeillll998tl . 

Corresponding variations are expected also for the 
CR pressure, P cr = | c f Q (3{p) g(x,p,t)pdp (/3 ~ 
1 for relativistic CRs, and (3 ~ p/(mc) for non- 
relativistic CRs) and for the CR ionization rate Ccr = 
C c f 0 v(p) g(x,p, t) (Ji (p) dp where a ;(p) is the cross sec¬ 
tion for the ionization of a hydrogen atom by a CR 
particle of momentum p, and the factor Cf accounts 
for heavy CR particles a nd for secondary electrons (eg 
(Soitzer fc Tomaskoj|!1968 1. These integrals will be com¬ 
puted elsewhere. Here we only stress that P CI and £ C r 
clearly increase with n cr . If they are nearly propor¬ 
tional to n cr , we then expect them to be nearly pro¬ 
portional to nf /2 as well, based on equation ©. The 
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values of P CI and £ cr should also drop when the density 
is large enough to cause the free-streaming of the CRs 
(when the second term on the r.h.s. of equation JS| is 
large). Such drop should occur at a gas density above 
« 500(i?o/10 fiG) cnr -3 in mostly neutral diffuse re¬ 
gions, and above ~ 10 5 (i?o /10 ^ G) cnr -3 in dark molec¬ 
ular cores, where rq is much smaller than in diffuse re¬ 
gions. The precise value of this critical density for free- 
streaming depends on the magnetic field strength and on 
the normalization of the CR spectrum. Furthermore, P cr 
and £ cr should decrease by a factor of 10-50 in the tran¬ 
sition from diffuse regions to dark molecular cores, also 
due to the reduced value of n\. 

Our result accounts for the heretofore unexplained en¬ 
hancement of the CR ionization rate in Q Persei diffuse 
gas (McCall et al. 2003; see also Le Petit et al. 2004) 
compared to dense molecular clouds, where i t is esti¬ 
mated using molecular abundance ratios l| William s et alJ 
flfiOSt IDotv et Ml 12002;: Pado an et al] 12004. _and refer¬ 
ences therein) or from the HI/H 2 ratio ^Goldsmit h Hr. Til 
120051 see Liszt (2003) for a discussion of other evidence 
concerning the CR ionization rate). Density estimates 
for the Persei gas using a variety of techniques are in 
the range 100-400 cm -3 , putting that gas near the up¬ 
per limit of our predicted CR density for diffuse gas, 
just before the free-streaming regime. In molecular re¬ 
gions, shielding from UV radiation allows carbon to re¬ 
combine and exist in neutral form or in molecules, result¬ 
ing in much smaller ion density then in diffuse regions, 
despite the large total gas density. This ion density is 
low enough that self-confinement is effective (up to a 
density of approximately 10 5 cm -3 ), but with a larger 
streaming speed (larger ion Alfven speed) than in dif¬ 
fuse regions. As a result, the CR density can be 10-50 
times lower than in diffuse regions, due to the constant- 
flux constraint expressed by equation (see Figure 1). 
Figure 1 should not be interpreted as predicting a one- 
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to-one or universal relation. For example, we expect a 
large dispersion in the relation between magnetic field 
strength and gas density and the CR flux normalization 
may vary with position relative to nearby CR sources. 

We have not discussed ionization due to electrons 
because the electron spectrum is very uncertain (see 
ICasadei fc Bindl2004ll . especially below 100 MeV, where 
it is sensitive to the models us ed to demodulate the CR 
electron flux at the Earth (e.g. IWebbeilll99 81 or used to 
disentangle the 7 -ray bremsstrahlung, inverse compton 
and unresol ved point source emission at MeV energies 
(e.g. lStron5ll2001 h However, CR electrons, like CR pro¬ 
tons, arc confined to magnet ic waves and should there - 
fore follow the CR protons llMelrose fc Wentzel 11970 1. 
and have the same density variations, although the crit¬ 
ical densities for free-streaming may be different. 

This work points out the possibility of significant spa¬ 
tial variations in the CR pressure, ionization and heating 
rates. However, here we have computed explicitly only 
variations of n cr , and not of P CI and £ cr . An explicit 
derivation of the corresponding variations in P cr and () cr , 
the inclusion of electrons, and a detailed discussion of the 
observed and predicted CR ionization rates will be given 
elsewhere. The variations we find should also be impor¬ 
tant for the ionization fraction in planetary atmospheres 
and protoplanetary disk ionization and chemistry, as well 
as the heating rate and CR -initiated ion-molecule chem¬ 
istry in molecular clouds. 
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